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We apply perturbative renormalization group theory to the symmetric phase of a dilute interacting 
Bose gas which is trapped in a three-dimensional harmonic potential. Using Wilsonian energy-shell 
renormalization and the e-expansion, we derive the flow equations for the system. We relate these 
equations to the flow for the homogeneous Bose gas. In the thermodynamic limit, we apply our 
results to study the transition temperature as a function of the scattering length. Our results 
compare well to previous studies of the problem. 
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O ■ I. INTRODUCTION 
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Dilute, repulsively interacting Bose gases in a three-dimensional harmonic trap have been at the focus of attention 
. of the theoretical physics community since the experimental realization of Bose-Einstein condensation (BEC) in 1995 
[H • A vast amount of literature has been published recently on this topic based on the mean-field approach. 
It is however known that mean-field theory is not appropriate near the critical temperature because the fluctuations 
dominate the mean field in this region. This follows from applying the so-called Ginzburg criterion 0, [f| and has 
triggered several attempts to go beyond mean-field theory near the critical regio n by means of various renormalization 
group (RG) techniques in the homogeneous interacting Bose gas 0, 0, H, 0, El Q ■ This recent work was added 
to a significant amount of existing literature on the application of RG methods to interac ting Bose gases, that was 
written when such systems were only an interesting theoretical problem 0, 0, [TEl ll6L Il8| . 

All these studies so far concern the homogeneous case; the real challenge however is to apply renormalization 
techniques to the experimentally relevant system, the trapped Bose gas This challenge is met here for the first 
' time, at least when we approach the transition region from above the critical temperature (symmetric phase). In 
this case the theoretical treatment is simpler than in the symmetry-broken phase and allows us to see how to apply 
renormalization techniques. In particular, we shall use a method similar to what in the homogeneous case is known 
as momentum-shell renormalization 0, [g, Il9| . In the trapped case it is energy shells instead of momentum shells 
that we will be integrating out. Of course, since we are dealing here with a finite system, there is, strictly speaking, 
no phase transition However, as we approach the thermodynamic limit, a quasi-phase transition develops, and 
renormalization methods are expected to enable us to study non-universal properties such as the critical temperature 
q and its dependence on the scattering length. 

The paper is organized as follows. Section II develops the theoretical methods necessary to apply the renormalization 
procedure to the weakly interacting, trapped Bose gas. Starting from a high-energy cutoff, we successively integrate 
out energy shells thus perturbatively creating an effective action for the low-energy field. This effective action is then 
cast into the form of the original action. In Sec. Ill, we derive the renormalization group equations for the chemical 
potential, the interaction, and the grand-canonical thermodynamic potential. We also refine our results by means of a 
first-order e-expansion. It is shown that, in the thermodynamic limit, the flow equations for the chemical interaction 
and the interaction have the same form as for a homogeneous interacting Bose gas. All nonuniversal properties, 
that are due to the presence of the harmonic trap, enter through the flow equation for the thermodynamic potential. 
Section IV is devoted to the study of a particular non-universal property, namely the transition temperature and its 
dependence on the scattering length in the thermodynamic limit. After a short summary in Sec. V, we discuss in 
Appendix A how the thermodynamic limit affects the relevant density of states entering the renormalization group 
equations of Sec. III. Appendix B contains a discussion of numerical methods we used for obtaining the results 
presented in Sec. IV. 



II. EFFECTIVE ACTION FOR THE LOW-ENERGY FIELD 

The grand partition function for the interacting dilute Bose gas can be expressed as a functional integral |2(], , 



i.e., 



z= 8[4>,<t>*] e -* s l+'+ m l (1) 
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We work in a /^-dimensional space, where the bosonic fields 4>(t, x) and <fi* (r, x) depend on spatial coordinates 
x = Xjj) and on imaginary time t. The fields are periodic in r with period h/3, where (3 = is the 

inverse temperature and ks denotes Boltzmann's constant. Assuming s-wave repulsive interactions, the Euclidean 
action characterizing the functional integral of Eq. Q is 



hp 



dr / d D x\cf)*(T,iz) 



(r,x) + ||0(r,x)| 4 
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where fj, is the chemical potential, m the particle mass, and the spatial coordinates are integrated over all space. The 
coupling constant for the short-range, repulsive interaction potential is denoted by g. For a dilute Bos e gas, i t has to 
renormalize to the two-body T-matrix when all two-body scattering processes are taken into account p. l2lll22| . Thus, 
in three spatial dimensions, at zero energy and in lowest order perturbation theory, this coupling constant is related 
to the s-wave scattering length a by g = Anah 2 jm. An improved relation for the scattering length a is discussed later 
[compare with Eq. (|26|l ] . 

We assume that the interacting Bose gas is trapped in a D-dimensional isotropic harmonic oscillator potential 
with frequency oj, i.e., V(x) = ^muj 2 x 2 . Furthermore let us define the trapping potential along each of its spatial 
dimensions as Vj(xj) = mLo 2 x 2 /2, j = 1, . . . , D. The corresponding normalized eigenfunctions are denoted by ip nj (xj), 
rij > 0, and have eigenenergies E rij + Eq/D, where E n . — hurij and Eq = Dhiu/2. The Z)-dimensional eigenfunctions 
are given by ipm,...,n D ( x ) = [x\) ■ ■ ■ ip nn (%d) with eigenenergies E n + Eq = hun + Eq where n = n\ + . . . + njj. 

We now wish to apply a modification of the momentum-shell RG procedure to the trapped interacting Bose gas 
characterized by the effective action of Eq . (J2J. In the homogeneous case the momentum-shell method has been 
explored extensively 0, 0, H H EHH E3 . In the case we are dealing with here, the noninteracting Hamiltonian 
(g = 0) derived from J3J does not commute with the momentum operator because of the presence of the trapping 
potential. We will therefore use harmonic oscillator eigenstates instead of momentum eigenstates and integrate out 
small energy shells. 

For the purpose of evaluating the partition function Q), we impose a high-energy cutoff by assuming that the 
maximum value that n can take is some large integer % > 1. We define an energy shell as the shell between — 8n^ 
and nA, where Sn\ is an integer such that 5n\/n\ < 1. In addition, we split the bosonic field into low-energy and 
high-energy components denoted by (j> < and </>>, respectively, so that 4>(t, x) = 0<(r, x) + 0>(t, x). Thus, the total 
bosonic field expanded on eigenfunctions of the noninteracting Hamiltonian is 
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where p l a — 2irl/Ti(3 are the Matsubara frequencies, <j>(l, n\, . . . , no) are the complex-valued expansion coefficients, and 
the prime on the sum over rii, . . . , no signifies the constraint n± > 0, . . . , no > 0, n± + . . . +Ud = n. Correspondingly, 
the low-energy field is 



oo riA— Sua 



E E E 
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(x) <p(l,ni, . . . ,n D ), 



and the high-energy field is 
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(x) ^(Z,ni, ...,n D ). 



The first stage of the RG procedure is usually referred to as Kadanoff transformation and consists of two steps; 
first, an effective action for the low-energy field is derived, and subsequently this effective action is cast in the form of 
the original action. For this purpose we proceed to the one-loop calculation of the effective theory for the low-energy 
field by integrating out the high-energy field. Analogously to the homogeneous case 13j, we thus obtain an effective 
action of the form 



Se ff [0<,0*<] = S[0<,0*<] + ^Tr[ln(l - G>£) 



(3) 



The "Tr" symbol denotes the trace in both the functional and the internal space of G'*!!. These latter quantities 
denote the Green's function for the high-energy field, i.e., 



G>(p ,H ho ) = 



B{p Q ,H ho ) 

B*(p ,H ho ) 
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and the self-energy for the low-energy field, i.e. 



1 ' j "H ^<(f,x)^(f,x) 20*(f,x)0 < (f,x) 



B(p ,H ho )= : 1 , ff ho = |- + V(x)- M . (4) 



with 



The hat indicates that these quantities are Schwinger-Fock operators |23l I24L l25| . 

As in the homogeneous case, we now perform an expansion of 10 up to second order in G > S. This expansion is 
particularly well justified in the context of the RG procedure used here, because the self-energy is multiplied by the 
high-energy propagator which is inversely proportional to the largest energy scale of the system, the cutoff energy 
fiumA. Thus an RG perturbative expansion is expected to have a wider range of validity than ordinary perturbation 
as employed, for example, in the mean-field theory context. Furthermore, the truncation at second order, that is at 
quartic interactions, is self-consistent with the truncation at quartic interactions of the original action (J2J. Higher 
order terms in powers of the low-energy field are discarded exactly as in the homogeneous case. A further discussion 
of these terms that we neglect here can be found in [3, 0] ■ The second-order expansion yields 

Tr[ln(l - G>±)} « Tr[-G>£ - ±(G>E) 2 ]. (5) 
Performing the sums over the Matsubara frequencies we find for the first trace 

-ft/3 r "A 
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and for the second trace 



/.ft/3 c- ™ A _L, 

Tr[G>tG > t}= I dr J d D x |^<(r,x)| 4 2g 2 £ £ f 2 (E n + E )\^ ni ... nn (x)| 2 , (7) 



where 



ME) = l + 2N BE (E), 

h(E) = ±(3N BE {E)[l + N BE (E)} + 1 + ™ BEi ? ) , (8) 

2{E-n) 

and N BE (E) = [e^ E ~^ - l]" 1 is the Bose-Einstein distribution. 

We proceed to the second step of the Kadanoff transformation noting that the first trace, Eq. JJJJ, is quadratic and 
the second trace, Eq. (JJJ, is quartic in the modulus of the low-energy field. Therefore the first trace can be interpreted 
as a correction to the quadratic part of the original action (|2j) , 

"A 

d( t i-V) = -g Yl E /i(£n + £o)|^„ 1 ...n D (x)| 2 , (9) 
and the second trace as a correction to the quartic part of (0), 

TIA 

dg = -g 2 ]T E f2(E n + E )\ip ni ... nD (x)\ 2 . (10) 

We observe that the corrections (jHJl and (|10J) are x-dependent and therefore not of the form of the original action 
@. However, in our subsequent treatment we are mainly interested in the limit of small trapping frequencies, i.e., 

[3huj « 1. In this case, the thermal wave length X t h = \J 2i:fi 2 (3 /m is small in comparison with the characteristic 

extension L = y/h/muj of the ground state of the harmonic trap. Therefore, it is expected that the thermodynamic 




FIG. 1: Schematic for integrating out a high-energy shell. Shown is a one-dimensional cut through the potential along the 
coordinate axis Xj. The high-energy eigenfunctions are located around Ea 3> -Eo- In Eq. 11 1 1 ^ - the trap potential Vj is thus 
negligible compared to Ea around the trap center. 



behavior is dominated by the properties of the interacting Bose gas in the center of the trap and boundary effects 
are negligible. In addition, this limit is consistent with the thermodynamic limit which we will concentrate on later 
and where ui and the particle number N tend to and oo, respectively, such that Nto 3 remains constant. Thus, let 
us focus on the region near the center of the trap which means Vj(xj) <C E n . + Eq/D, j = 1, D. FigureHshows a 
schematic of the situation under consideration. We now employ the JWKB approximation for the harmonic oscillator 
eigenfunctions [1^,|23, i.e., 



^n, (Xj) 



^(2m) 



1/4 



[E^.+Eo/D-Vjixj)] 



1/4 



cos 



dy (j2m[E ni + E /D - V 3 (y)] - J 



(11) 



with = Ly / 2nj + 1 and = —Ly/2rij + 1 denoting the right and left classical turning points. Expression (|1 If) 
is valid in the classically allowed region < Xj < except for a small area near the turning points. In particular, 
near the trap center it is possible to neglect the trapping potential Vj(xj) in the denominator of (|llfl . Furthermore, 
in this region the high-energy eigenstates oscillate very rapidly on the length scale set by L. Therefore, we can safely 
approximate each of the squares of the cosines appearing in 10 and (|10|) by their spatial averages 1/2. The corrections 
in Eqs. © and J7J now assume the forms 



™ A (9mtA D / 2 

n—n/y — 8n\ \ / 



(12) 



and 
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(13) 
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where 

£ / = ( 14 ) 

m...n D w (ni + 5 )...(n.D + 5) 

We note that because we are focusing on the center of the trap, there is no explicit dependence on the trapping 
potential anymore on the right-hand side of i|12fl and 113(1 . We can therefore interpret the right-hand side of l|12[) as 
a correction to the chemical potential only and set dV = 0. In other words, the trapping potential and consequently 
the trapping frequency receive no non-trivial corrections and exhibit trivial scaling only. 



III. RENORMALIZATION GROUP EQUATIONS 
A. Derivation of the flow equations 

At this stage we replace the sum over n in Eqs. I|12|) and (|13fl by an integral over n. It is well known that, as long as 
the functions we are summing over are smooth, this is a satisfactory approximation |26l l28t l29 | . and it is particularly 
good in the regime we are considering here, (3fioj <C 1 |30|. We can now consider the renormalization step to be 
infinitesimal and pursue the analogy with the homogeneous calculation. Therefore, we shift the origin of integration 
over n in (|12|l and l|13[l . i-e. n — > n — -j, thus eliminating the zero-point energy from the functions /1, fi in JSJl. 

In order to emphasize the analogy to the case of a homogeneous interacting Bose gas we introduce a generalized 
momentum vector p = (pi, ...,pjj) through the relation 

P 2 

nuj nj = ^, Pj >0, j = l,...,D. (15) 

Obviously, the modulus p of this generalized momentum is given by hum = p 2 /2m. We can also define a cutoff hA 
for the generalized momentum vector by the relation hu>n A — (hA) 2 /2m — E A . Correspondingly, 1/A may be viewed 
as the smallest length in our problem, i.e., A~ x <§; \ t h or equivalently /3a <C /3 with (3 A = m/h 2 A 2 — 1/2E A . As in 
the homogeneous case, we may parametrize the generalized momentum in terms of this cutoff and of a dimensionless 
continuous parameter / according to p/h = Ae~ l . Thus, differentiating Eqs. I(12|l and (|13|) with respect to / we obtain 
the flow equations 

^ = -g (A e- l ) D d(l) h(E A e- 21 ), 

f = -g 2 (A e- l ) D d(l)f 2 (E A e- 21 ) (16) 

with 

d(l) = ^d(LA e- l f- D p[{LA e~ 1 ) 2 /2 - D/2}. (17) 

This latter factor describes all effects originating in the discrete energy level structure in the isotropic harmonic trap. 
In the case of a homogeneous interacting Bose gas it assumes the constant value dh = 1/2-k 2 for the physically relevant 
case of D = 3. 

The effective action of Eq. Q can now be cast in the form of the original action. This is achieved by reinstating 
the cutoff for n to its original value n A , or, equivalently, the cutoff of the modulus of the generalized momentum q 
of the low-energy field to its original value hA (note that we are using q for the modulus of the momentum of the 
low-energy field as opposed to p which is the modulus of the momentum of the high-energy field). To this end, we 
rescale each component qj of the momentum vector q = (qi, qo) of the low-energy field according to qj(l) = qje 1 . 
Thus, the rescaled modulus of the momentum vector of the low-energy field reads q(l) — qe l and it is cut off at HA. 
This trivial rescaling together with the demand that the effective action (01 remains formally the same after each 
renormalization step induces the trivial scaling of the parameters of the effective action, x(V) = xe~ l , t(1) — re~ 21 , 
fi(l) — fie 2l ,g(l) — ge^ 2 ~ D ^ 1 , (311) = (3e~ 21 , </>(/) = (f>e Dl ^ 2 . We observe that all quantities scale as in the homogeneous 
case. In addition, we define here the trivial scaling of the oscillator frequency, ui(l) — uie 21 . Consequently, the trivial 
scaling of the oscillator length is given by L(l) = Le~ l . It is convenient to introduce the dimensionless parameters 
b{l) = (3(l)/f3 A ,M(l) = p Af i(l),G(l) = p A A D g(l)/b(l),Q(l) = Tl(3 a uj{1) = (LA e~ l )~ 2 and the dimensionless cutoff 
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energy E > = (3\E\ — 1/2. We want to point out that the dimensionless quantities b(l) and M(l) scale trivially as 
the corresponding dimensional parameters (3(1) and fi(l), whereas Gil) scales trivially as G(l) = Ge 1 - ^ 4 ^ 1 . In terms of 
these quantities the renormalization group equations for the dimensionless chemical potential and the dimensionless 
coupling strength are finally given by 

= 2M(l) - G(l)d(l)b(l) [1 + 2N(l)] , 

^ = (4 - D)G(l) - G(lf d(l)b(l) {4 6(0 N(l) [1 + N(l)} + 2( ^^ )} } (18) 

with the scaled Bose distribution 

N(l) = {e b ^ E >- M ^-iy 1 . (19) 

We observe that there are two differences from the homogeneous RG equations. 

The first difference is that the maximum number of RG steps in the trapped case is finite. This is due to the presence 
of the zero-point energy which has no homogeneous counterpart. The number of renormalization steps I required to 
reach an energy scale E{p) = p 2 /2m = hum is defined by E(p)e 21 = E\ and therefore I = In y/ E\/ E(p). This means 
that, if -Emm is the minimum value the energy can take, the maximum value of RG steps is I* = y E\/E m i n . In the 
homogeneous case E m - m = and I* = 00 whereas in the trapped case E min = Dhto/2, which results in a finite number 
of steps, 

I* = In - 1 = In - 1 (20) 

Of course, this maximum number of steps must be consistent with the requirement that throughout the renormalization 
procedure the Bose-Einstein distribution of Eq. (|19|l must remain positive so that the particle density remains positive. 

The second difference is due to the fact that in the trapped case the energy spectrum is discrete. Although we 
replaced the sum over n by an integral over n we kept the sums in the expression for p(n) of Eq. (|14|) and thus retained 
part of the discrete nature of the problem. This method is also employed in treatments of the non-interacting trapped 
Bose gas [2^, |2^ where it is also shown that, for D = 3, the corrections due to the discreteness play a significant 
role only for small particle numbers, typically N < 1000. The consequence of retaining the sums in p(n) is that the 
quantity d(l) of Eq. lfT7| appearing in the RG equations (fl%|l differs from the homogeneous case in two respects. First, 
it depends on the renormalization step I whereas the corresponding homogeneous quantity is a constant, e.g., for 
D = 3 we have d^ = 1/27T 2 . Furthermore, d(l) does not coincide with the density of trapped states, whereas d^ does 
coincide with the density of homogeneous states 0, 0] . 



B. Thermodynamic limit 

In our subsequent discussion we are particularly interested in the thermodynamic limit of the RG flow equations. 
As mentioned above, in the case of an isotropic harmonic trap this limit is defined by letting the trapping frequency 
to tend to zero and the number of particles N to infinity in such a way that Ncu 3 remains finite. It is shown in 
Appendix A that in this limit the quantity d(l) approaches the constant value 1/27T 2 in three spatial dimensions. We 
thus conclude that in the thermodynamic limit the RG flow equations l|18H for the trapped interacting gas assume 
the same form as in the homogeneous case. This is one of the main results of this work. In particular, it follows that 
the universal critical properties are not influenced by the presence of the trap. As we will see below, all effects on the 
(nonuniversal) thermodynamic properties that originate from the isotropic harmonic trap are contained in the flow 
equation for the grand thermodynamic potential of Eq. Ij24j) below. It involves the density of states in the trap and 
therefore differs from the corresponding flow equation for a homogeneous Bose gas. 

It is also instructive to study the transition to the thermodynamic limit in terms of the flow resulting from the RG 
equations l)18[l. Figs.GIa) and (b) compare typical flow trajectories close to and further away from the thermodynamic 
limit. The proximity to this limit is determined by the value of f2(0) = V(-^A) 2 which differs by several orders of 
magnitude between the two plots. Within each of the diagrams, however, the trajectories were obtained for fixed 
initial values of b, G, and fl, but differing values of the scaled chemical potential M(0). The dashed curves show 
representative trajectories for the corresponding homogeneous flow, i.e, O(0) = 0. These trajectories clearly reflect the 
presence of the unstable fixed point in the homogeneous flow which is indicative of a second-order phase transition. 
The trapped flow of Fig. [3a), which has a very small but finite value of 0(0), does not have a fixed point anymore in a 
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FIG. 2: Flow trajectories determined from the trapped RG equations 1181 1 (full curves) and, in the thermodynamic limit 
(lo — > 0, N — *• oo, Nlo 3 = constant), the corresponding homogeneous trajectories (dashed curves). All trajectories have the 
same initial values of 6(0) = 10657 and G(0) = 4ir but different values of M(0). For the homogeneous trajectories Af(0) varies 
between 2.1387 x 10" 6 and 2.1391 x 10" 6 . The corresponding value for aN 1/6 /L is 6.12 x 10" 3 . The trapped trajectories in 
(a) are for Q(0) = 1/(LA) 2 = 10~ 9 and M(0) between 2.1195 x 10~ 6 and 2.128 x 10~ 6 . The trapped trajectories in (b) are for 
fi(0) = 10~ 6 and M(0) between 1.4 x 10" 6 and 2.1 x 10~ 6 . The quantity G is related to G by G = Gd(l). 



strict sense, but it stills resembles very closely the homogeneous flow. In this way, it is representative of a quasi-phasc 
transition, that is expected for systems close to the thermodynamic limit. In Fig. Hfb), fi(0) is larger by a factor 
of 10 3 , and the flow now differs significantly from the homogeneous one, although some characteristic features are 
still present. One should also note that M(0) now has to be changed over a much wider range than in (a) in order 
to produce the characteristic change in the large-^ behavior of the flow trajectories, i.e., switching from negative to 
positive final values M(l*). This is in accordance with the expectation that the phase transition becomes more and 
more smeared out, as one moves away from the thermodynamic limit. 

We now turn to the calculation of the grand thermodynamic potential that allows to calculate all thermodynamic 
properties of the trapped gas. In analogy with the symmetric phase of the homogeneous gas, within the quadratic 
approximation [3l| , the grand thermodynamic potential in the symmetric phase of the trapped gas is given by 



w = - Y In [l - e-KSo+Bo-ri] =-Y T In 

hi,...,jid n=0rii...nD 



_ -p(E n +E -fJ.) 



(21) 



In D — 3, the constrained sum in the above equation can be determined easily 28, 29, 30] and w reduces to 



w 



ft 
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n=0 L 



n 3n 



In 



-p{E n +E a -lt) 



(22) 



Replacing the sum over n with an integral in (|22|) , shifting the origin so that the zero-point energy is eliminated from 
the integrand, and using the generalized momentum p — h,Ae~ l instead of n (compare with Sec. Ill A), we obtain in 
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the thermodynamic limit 

1 /' x „ I 



w 



dl' „ In 

l'=0 



p J v=0 



e -b{l')(E > -M{l')) 



(23) 



It is convenient to define an 'intensive' dimensionless thermodynamic potential W(l) — H(0) 3 Paw(1) for which the 
scaling of the extensive potential w in the thermodynamic limit is appropriately taken into account. Differentiating 
with respect to the continuous variable I we finally obtain the flow equation 

We note that in the above RG equation the chemical potential depends on I and therefore, through (JTSJ, on the 
interaction. In other words, the interaction enters the RG equation for the grand thermodynamic potential implicitly 
through the non-trivial scaling of the chemical potential, which is the main difference to the case of a trapped ideal 
gas. In order to determine the grand thermodynamic potential, we have to solve Eqs. (|18(l and 1|24|1 with initial 
conditions W(0) = 0, M(0) < 1, G(0) = G(0)6(0) > 1, b(0) > 1, and fi(0) < 1. 

The renormalization group equations (|18ll can be refined in a simple way by means of a first order e-expansion 
|l9l l32] | . This method is expected to yield better quantitative results in particular as far as critical properties are 
concerned. For this purpose we set 4 — D = e in (|18|l and formally treat this quantity as a small parameter. We then 
expand systematically all quantities on the right hand sides of the RG equations l|18fl around (M, G) = (0, 0) up to 
second order in this small parameter. Thereby one assumes that the quantities M and G remain of the order of e. 
Thus, we expand the first equation of l|18|) up to first order in M and the second equation of l|18|l up to zeroth order 
in M. As a consequence, in the thermodynamic limit we find the modified RG equations 

= 2M(l) - G(0t^Td K0 {2A m [6(Z)£>] + 2b(l)N m [b(l)E > ](l + Ar m [&(Z)£>])M(Z)} , 

^ eG{l)-G{lf^b[l)x 



dl w w (2tt) j 

l + 2N m [b{l)E > ] 



46(0 AU &(/)£>](! + AUK0£>D \ (25) 



2E> 

with Dd denoting the surface of a Z?-dimensional unit sphere (compare with Appendix A) and with the modified 
Bose-distribution N m (x) = Nbe(x; fi = 0) = [exp(x) — 1] . In the framework of the e-expansion one sets D = 4 on 
the right hand side of Eqs. I|25|l . In three spatial dimensions these RG equations have to be solved for e — 1 together 
with the flow equation (|24|l of the grand thermodynamic potential. 



IV. CRITICAL TEMPERATURE 



Solving Eqs. Ijl8|l or Q25JI together with Eq. I|24(l allows us to determine any thermodynamic property of an in- 
teracting, trapped Bose gas at and above the BEC phase transition. First of all, let us summarize the results for 
the critical exponents which have been derived elsewhere |6|, l32t |33j in connection with the study of the homoge- 
neous flow. In the thermodynamic limit, the unstable fixed point of the flow equations (|18|l with D = 3 is given by 
(M*,G*) = (1/12, 57r 2 /72). The critical exponent for the correlation length is v = 0.532. Using first order e-expansion, 
this fixed point is shifted to the position (M*,G*) = (e/10, e7r 2 /10) and its eigenvalues are given by Ai = 2 — 2e/5 
and A2 = — e with e = D — 4 = 1. As a consequence the corresponding critical exponent for the correlation length 
assumes the value v — 1/Ai = 0.600 which is significantly closer to the experimental value of v = 0.670 [33| . 

In the following, however, we wish to focus our attention on the dependence of the critical temperature on the 
interaction strength in the thermodynamic limit. The numerical investigation of critical properties requires us to find 
the critical flow trajectories in the (M, G, b) space, i.e., those trajectories that asymptotically reach the fixed point for 
I — > 00. In Appendix B, we briefly comment on strategies for accomplishing this task. There exists a continuous family 
of different critical trajectories, the so-called critical manifold, and each trajectory can be assigned a specific value 
for aA 1 ^ 6 /L, i.e., a suitably normalized scattering length, as we see below. Scanning the different critical trajectories 
thus allows us to explore a wide range of scattering lengths. The thermodynamic potential W at criticality is then 
determined by integrating Eq. (|24|l along a critical trajectory. Thereby, the initial and final values [M(0), G(0), 6(0)] 
and [M(lf),G(lf),b(lf)} of the trajectory need to be pushed sufficiently outwards so that the calculation of W is 
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converged. Typically, this requirement implies that one has to choose M(0) <C 1 and G(0)6(0), b(0) 1. In this way, 
it is ensured that the cutoff A corresponds to the largest energy in the problem, and relevant physical quantities such 
as aiV 1 / 6 /L and the critical temperature become cutoff- independent. 

Applying this procedure, it turns out, however, that for the largest scattering lengths shown in the figures below, 
convergence of W cannot be achieved even if the initial conditions are pushed very far outwards. The results thus 
become somewhat cutoff-dependent. To explain this behavior, we argue that in this regime our RG model effectively 
describes a gas of hard spheres. We then attribute the cutoff-dependence to the fact that for strong interactions, 
the critical temperature no longer depends only on the scattering length, but is also sensitive to finer details of 
the potential. To establish the connection to a hard-sphere gas, we note that the bare two-body scattering length 
pertaining to a set of given initial conditions is obtained from the flow equation for g(l) = g{l)e^ D ~ 2 ^ 1 in the limit of 
zero temperature and zero chemical potential. The quantity g(l) characterizes the renormalized interaction strength, 
after the trivial scaling has been removed. Identifying the s-wave scattering length a with the renormalized interaction 
strength g(oo) according to g{oo) — Aixah 1 jra we obtain the relation 

aA= 6 y (26) 

for D = 3. This latter relation for the zero-energy s-wave scattering length is in agreement with the Lippmann- 
Schwinger analysis for the T- matrix Q. We then compare this result to zero-energy two-body scattering with 
an interaction potential V^(x — x') = <3(R — |x — x'|)?i 2 KQ/2m r , i.e., a finite-height repulsive potential of radius R 
[35j . Thereby m r = m/2 denotes the reduced mass and kq is a measure of the potential strength. If n R ^> 1, the 
corresponding scattering length a is given by a/R— k R/(1 + k R). With the help of the identifications R = ir/2A 
and KoR = G(Q)b(0) /2tt 2 , this expression can be mapped onto Eq. (|26(l . In this way, we can relate the interaction 
used in @ to scattering between quasi-hard spheres, and 7r/2A is interpreted as the sphere radius. 

As emphasized in Sec. Ill B, in the thermodynamic limit the flow equations (|18f) and I125H reduce to the equations 
of an interacting homogeneous gas. Thus, all effects originating from the presence of the isotropic harmonic trap 
are contained in the flow equation (|24Jl for the dimensionless 'intensive' thermodynamic potential W(l). As W(l) 
only depends on M(l), but not on G(l), the particle interactions affect the thermodynamic properties only through 
the non-trivial scaling of the dimensionless chemical potential M(l). In the special case of negligible interaction, 
i.e. G(l) = 0, M(l) scales trivially, i.e., M(l) = M(0)e 21 , so that Eqs. (0 and ||2U) reproduce all thermodynamic 
properties of a non-interacting Bose gas in an isotropic harmonic trap within the framework of a grand canonical 
ensemble. 

For a given critical trajectory with initial conditions [M(0), G(0), 6(0)], the calculation of the corresponding crit- 
ical temperature T c proceeds as follows. With the critical trajectory, we associate the scaled particle number 
s = —dW/dM\ b ^ = n 3 N. A way to accurately calculate this quantity numerically is outlined in Appendix B. 

An ideal gas with the same value of s = Q 3 N has a critical temperature Tj? that is determined by P\k B T^ = [sC(S)] 1 / 3 
with ((x) the Riemann zeta function [3(| . The ratio T c /T® can thus be expressed as 

T c /T c ° = C(3) 1/3 /(s 1/3 6). (27) 

In our subsequent discussion we want to investigate the dependence of the critical temperature on the s-wave scattering 
length of the interacting Bose gas. From Eq. I|26|l we find that a particularly convenient measure for the scattering 
length is the dimensionless parameter 

aN^/L = s^aA = ,1/6 G(0)b(0) 

4^+|G(0)6(0) 

For each critical trajectory, we calculate T c /T® and aN 1 ^ / L, and in this way we obtain the dependence of the critical 
temperature on the scattering length. 

Figure [3] summarizes our main numerical results for the critical temperature. Figure EUa) shows the RG calcula- 
tions according to Eqs. lf2"5t (bold curve) and ifTSI) (dashed) together with the mean field (MF) approximation j3|J 
(dotted) and the predictions of Rcfs. [23] (dot-dashed) and £3] (dot-dot-dashed). It is remarkable that all results 
are within reasonable quantitative agreement. In particular, both RG results, which match rather closely for the 
values of aiV 1 / 6 /L shown in Fig.|^a), lie above the MF approximation. This is generally expected [38| as the critical 
fluctuations, that are neglected in the MF theory, should lead to an increase in T c . In the limit of aN 1/6 /L -> 0, the 
RG curves tend to approach the MF result, and the ratio (T* G - T C MF )/(T C ° - T C MF ) seems to assume a constant 
value of about 7%. The prediction of [3^] closely agrees with the RG curves, whereas the result of is intermediate 
between MF and RG. 
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FIG. 3: Critical temperature as a function of the scattering length for the harmonically trapped interacting Bose gas. Bold 
curve: RG result after e-expansion according to Eqs. 11251 . dashed: RG result without e-expansion according to E qs. 11181 . 
dotted: mean field approximation [3^1 . The dot-dashed and dot-dot-dashed curves in (a) show the predictions of Refs. IS?! and 
[3S|l . respectively. 



In Fig.|3Jb), we have extended our RG calculations to larger values of aiV 1 / 6 / L. As mentioned above, for aiV 1 / 6 / L 
larger than, approximately, 0.3, the quantitive results become somewhat cutoff-dependent. The general behavior of 
the curves, however, does not change. We now observe a profound difference between the treatments without and 
with e-expansion. The former predicts an increase in T c /T® for larger scattering lengths, whereas the latter indicates 
a decrease. On physical reasons, the second behavior appears to be more plausible, as for large particle interaction, 
an inhibition of Bose-Einstein condensation is expected. On these grounds, the e expansion seems to be more reliable, 
but more work is certainly required to support these results further. 

For the sake of comparison, we also show, in Fig.^ the results for the RG treatment of the homogeneous interacting 
Bose gas. In this case, the relevant thermodynamic potential Whom is obtained from integrating j|| [tj 



dWhom 1 



dl 2tt 2 6(0) 



-«l I1 [i_e- 6 ( | )^>-M(0]] ) (29) 



whereas the flow equations for M(l) and G(l) are still given by the thermodynamic limit of Eqs. I)25|) and Eqs. 
(|18fl . respectively. In Fig. 0] the scattering length is scaled to the cubic root of the particle density n p , so that 

aiip^ 3 = s^ m G(0)b{0) I '{At: + 2G(0)6(0)/tt) with s hom — —dWh om /dM\ b ^ Q, Q y The scaled critical temperature is 
obtained from T c /T° = [£(3/2)/sh om ] 2 / 3 /27r&(0). Together with the results for the symmetric phase we also display 
the calculation for the symmetry-broken phase as described in Ref. |j. For small an\j 3 it was established in that 
the computation in the symmetry-broken phase predicts a linear dependence of the critical temperature on a, i.e., 
T c /T c ° = Can],^ 3 with C ~ 3.4. For the symmetric phase without and with e-expansion, we also find an approximately 
linear behavior with constants C — 1.0 and 1.4, respectively. All these results fit reasonably well into the 'zoo' of 
current predictions |39| . In fact, the value for C obtained from the symmetric phase with first-order e-expansion 
matches rather closely the results of Refs. which may be regarded as the best current estimates |42^ . 



11 



1.3 




FIG. 4: Critical temperature as a function of the scattering length for the homogeneous interacting Bose gas. Bold curve: RG 
result after e-expansion according to Eqs. 1251 . dashed: RG result without e-expansion according to Eqs. 1181 . dotted: RG 
result for symmetry-broken phase [ji- 



lt is also of interest to consider the regime of large scattering lengths. In the symmetry-broken phase, the calculation 
predicts a monotonous increase of T c with the scattering length, a behavior which is probably unphysical. Remarkably, 
however, both curves for the symmetric phase show a maximum and a subsequent decrease in the critical temperature. 
The location and in particular the height of the maxima are somewhat cutoff-dependent, but altogether they do not 
differ too strongly from the results displayed in Ref . . As outlined above, we attribute the cutoff-dependence to the 
fact that in this regime our RG model effectively describes a gas of quasi-hard spheres. For strong enough interactions 
the critical temperature cannot be parametrized by the scattering length alone, but also depends on the finer details 
of the potential. Varying the cutoff, e.g., amounts to changing the radius of the hard spheres. 



V. SUMMARY AND CONCLUSIONS 

We have applied the Wilsonian momentum-shell renormalization group technique to the harmonically trapped 
interacting Bose gas. By integrating out small energy shells in the partition function and applying the e-expansion, 
we have obtained flow equations for the thermodynamic variables. In the thermodynamic limit, the flow equations 
for the chemical potential and the interaction potential reduce to the corresponding relations for the homogeneous 
interacting Bose gas. The presence of the trap becomes manifest only through the modified flow equation for the grand 
thermodynamic potential. From the flow equations, we have calculated the transition temperature as a function of the 
scattering and found our results in good agreement with previous approaches. In future work, we plan to extend the 
methods presented here to more general potential shapes in order to assess the general applicablity of our approach. 
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APPENDIX A 

In this Appendix, we discuss a method for calculating the constrained sum p(n) of Eq. I|14l) in the limit of large 
quantum numbers and wide traps, i.e., n» 1 and f3huj — > 0. In this case, the difference between the energy levels 
of the harmonic oscillator vanishes and we can therefore replace the sums in (I14|l by integrals. Furthermore, in the 
limit of large values of n the dominant contribution to Eq. (|14f> comes from values Ui, njj 3> 1 so that we obtain 
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FIG. 5: Numerical calculation of p(n) of Eq. 11141 (full curve) and approximations 2n^/n (dashed) and 2-K^/n — 5.701 (dotted). 



the approximate relation 

dn p(n) = dn dn\...dnr> 



Here we have converted the constrained integral over rij, j — 1, D and the integration over n into an unconstrained 
integration over rij and then changed variables to y/nj. Next we define the vector r = (yfn~i, y/n~o) which fullfils 
r 2 = ii\ + ••• + njj = n. This implies that 



where fljj = 2-n D l 2 /T(D/2) is the surface of a D-dimcnsional unit sphere, and r(x) denotes the Gamma function. 
The factor 1/2-° is due to the fact that the integration variables take only positive values. We have thus derived the 
main result 

il D n 2 
P( n ) = 2 

which is valid in the limit n — > oo. For D = 3 we have O3 = 47r and pin) — lnyfn^ so that in the limit £1(1) = 
1 / (LAe~ 1 ) 2 3> 1 the quantity d(l) of Ea. H17J) reduces to the result dh = 1/2tt 2 and the RG equations (|18[) reduce to 
the ones for a homogeneous interacting Bose gas in the absence of a trap. Figure |3] shows a diagram of p(n) together 
with the approximation It is also possible to compute numerically the next-order correction in the high-n 

expansion, 2-K^/n — 5.701, which is almost indiscernible from the exact result, except near n = (see inset). 



APPENDIX B 



In this Appendix, we briefly discuss some of our numerical methods. We first describe how we determine the 
derivative dW/dM \ b ^ which is required to calculate the quantities T c /T® and aA 1 / 6 /L of Eqs. l(2T)) and J2EJ), 
respectively. From these quantities, we obtain Figs. 01 and 01 showing the dependence of the critical temperature on 
the scattering length. 

A straightforward approach would consist in using finite differences, i.e., dW/dM \ b ^ q/ \^ [W(Mq + 5M) — 
W(Mq)]/8M with Mq the chemical potential for which the derivative is required and SM a small variation. However, 
this method turns out to not produce well-defined results, in particular for small values of aN 1 ^ 6 /L. Therefore, we 
apply a scheme which is closely related to the method of linear stability analysis frequently used in the theory of 
dynamical systems [44^. In this scheme, one propagates the flow trajectory, from which W(M ) is calculated, along 
with an infinitesimal linear deviation from the initial conditions. 
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To be specific, we consider the evolution of quantities SyX(l) with X, Y = M or G. These quantities are to be 
interpreted as 8yX(l) = 5X(l)/8Y(0), i.e., they describe the change SX(l) of a flow trajectory if the initial condition 
is changed infinitesimally by SY(0). If we write the flow equations l(T%)) or $ZJfy symbolically as 

"-f = Im(m, g), 

Jl = fc(M,G), 
then the deviations 5yX(l) obey the differential equations 




with the initial conditions 5m M(Q) = 1 and 5mG(0) = 0. The matrix of partial derivatives has to be evaluated 
along the original flow trajectory. Writing Eq. Ij24(l symbolically as dW(l)/dl — F(M), the required derivative 
dW/dM l^p-j (3( )= dW/dM g,, Q s (I — > oo) is finally obtained from 

d (dW \ _ dF 

dl [dM ) ~ dM dMM 

with the initial condition dW/dM |w % 0. We have verified that the results from this method agree with the 

finite differencing scheme in regimes where the latter is applicable, e.g., for larger values of aN 1 ' 6 / L. 

We now outline how to find the critical trajectories, i.e., those trajectories that asymptotically reach the fixed point. 
A systematic and reliable way to determine the critical M(0) for given [G(0), b(0)] consists in propagating trajectories 
with increasing M(0), starting from M(Q) = 0. For these trajectories M(l — > oo) — > — oo initially. However, for a large 
enough M(0) = Mi, one eventually finds M(l — » oo) — » 1/2. The crossover between these two opposite behaviors 
occurs at the critical trajectory. The critical M(0) is thus bracketed by and Mi and can now be found, e.g., via 
bisection. Another method, which is computationally more efficient but works well only for smaller scattering lengths, 
is based on backwards propagation from the non-trivial fixed point. To this end, one first finds the critical manifold 
of the fixed point in the linearization approximation. For a given b(lf) <C 1 and starting from points (M(Z), G(l)) on 
this critical manifold, very near the fixed point, we propagate the thermodynamic limit of Eqs. (|18fl or 125|) backwards 
up to a value I = U for which G{li) is sufficiently large so that it is ensured that aA < 1. In this way we find the 
initial conditions of a critical trajectory [M(7j), G(k), Further critical trajectories are computed by repeating 

the above procedure using different values of b(lf). 
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